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Abstract 

The purpose of this paper is to study the properties of kinetic models for traffic flow 
described by a Boltzmann-type approach and based on a continuous space of microscopic 
velocities. In our models, the particular structure of the collision kernel allows one to find 
the analytical expression of a class of steady-state distributions, which are characterized by 
being supported on a quantized space of microscopic speeds. The number of these velocities 
is determined by a physical parameter describing the typical acceleration of a vehicle and 
the uniqueness of this class of solutions is supported by numerical investigations. This shows 
that it is possible to have the full richness of a kinetic approach with the simplicity of a 
space of microscopic velocities characterized by a small number of modes. Moreover, the 
explicit expression of the asymptotic distribution paves the way to deriving new macroscopic 
equations using the closure provided by the kinetic model. 
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1 Introduction 

The purpose of kinetic theory for traffic flow is to provide an aggregate representation of the 
distribution of vehicles on the road, thanks to a detailed characterization of the microscopic 
interactions among the vehicles, which play an important role in the macroscopic trend of the 
flow. The goal is to obtain information on the macroscopic characteristics of the flow without 
assuming previous knowledge on the dependence of the mean velocity on the local density of 
traffic, as it is done in standard macroscopic traffic models. See in particular the prototype of 
macroscopic models, DU, or the reviews in mum. More refined macroscopic models consider 
a system of equations, instead of a single equation, see [2j, and/or they prescribe different flow 
conditions at certain stages, building phase transitions within the flow, 01 nans], but still it is 
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necessary to complete the model with a closure law, derived from heuristic or physical arguments, 
or from experimental data. Kinetic models provide quite naturally a closure law, which is linked 
in general to the equilibria of the kinetic model. Another way to derive closure laws is through 
microscopic “follow the leader” models HU2B]- For a review on the derivation of macroscopic 
traffic models from the microscopic “follow the leader” ones and from the mesoscopic kinetic 
theory, see [3]. 

Several kinetic approaches have been proposed, starting from the pioneering work of El 

and later m- These models were based on a Boltzmann-type collision term in which the 
cross section, giving the probability of an interaction between two particles, is replaced with 
a probability distribution depending on the local traffic conditions. The equation describes 
the relaxation of the kinetic distribution in time due to the acceleration and slowing down 
interactions among vehicles. However, the interaction integrals appearing in kinetic Boltzmann- 
type models for traffic flow based on a continuous velocity space, see EJ, typically do not provide 
the analytical expression of the equilibrium distribution and they are very demanding from a 
computational point of view. For this reason, two main approaches have been taken into account 
in order to compute the time-asymptotic distribution or to reduce the computational cost: on the 
one hand, one may consider Vlasov-Fokkcr-Planck type models in which the interaction integrals 
are replaced by differential operators, see uni m ns]; on the other hand, one may consider 
simplified kinetic models with a small number of velocities, namely the discrete-velocity models, 
see 7J0- 

In this work, following the classical Boltzmann-like setting of binary interactions, we study a 
kinetic model based on a continuous velocity space, which does not suffer from the aforementioned 
drawbacks. To this end, we focus only on spatially homogeneous problems and we investigate 
the interplay between the microscopic rules appearing in the collision kernel and the equilibrium 
solutions, finding analytical expressions of the time-asymptotic distribution, which result in a 
realistic macroscopic model for traffic flow. 

As in EH, our transition probability is characterized by the fact that drivers react to the 
presence of other vehicles, deciding whether to modify their speed according to the overall traffic 
conditions and to the particular velocity of the cars around them. Thus, the decision of whether 
and how to modify one’s speed depends on the local traffic density, or better, on the free space 
available, as we argue in [25] , Clearly, other choices were considered in literature. For example, in 
m the authors assume that drivers react to the local mean speed and they decide to accelerate 
or to brake by comparing their velocity to the speed of the flow. Here, the possible speeds 
available to the driver are naturally the driver’s current speed, and a set of speeds which depend 
on the velocity of the vehicles ahead and on the local density. In particular, the probability of 
accelerating increases monotonically with the free space available. Instead, the probability of 
braking increases as the road becomes congested. 

This framework permits to take the stochasticity of the drivers’ behavior into account, thanks 
to the probability distribution which assigns a weight to the possible driver’s decisions, while 
maintaining the general kinetic setting, based on the deterministic evolution of the distribution 
function. We propose two models that follow the framework just described: one is based on 
quantized velocity jumps, i.e., if acceleration occurs the new speed is obtained by increasing the 
pre-interaction velocity of a fixed quantity Aw (<5 model); the other one is based on a continuous 
uniform distribution defined on a bounded interval parametrized by Aw (% model), see also 
[ 11113 - 

In this paper, Aw is a finite parameter which models the physical velocity jump performed by 
vehicles when they increase their speeds as a result of an interaction. Clearly, this parameter may 
depend on the mechanical characteristics of vehicles, see Eli but in this paper we will assume 
that Aw is fixed. In the section on macroscopic properties, ^5] we show that Aw is related to the 
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Figure 1: Connection between the S model and a discrete-velocity model, having the same 
steady-state distribution. 


maximum acceleration and we discuss how this parameter can be chosen through experimental 
data used by Lebacque m- 

This paper also shows that the fundamental diagrams (or closure laws) obtained with the 
perhaps more natural, but more complex and more computationally demanding, \ model are 
very close to those provided by the simpler S model. We thus investigate the equilibria of the S 
model, both from an analytical and a numerical point of view. Analytically we find that velocity 
distributions formed by a linear combination of Dirac S ’s may be equilibria only if the S ’s are 
centered at velocities spaced by multiples of Av. Next, we compute equilibria using a numerical 
scheme capable of converging also to possible absolutely continuous equilibria. Here we find only 
the quantized equilibria described above, independently of the discretization parameters. This 
fact suggests that the class of discrete-velocity equilibria is the only one that the continuous- 
velocity Boltzmann-type S model possesses. This situation is summarized graphically in Figure]!] 

The paper is organized as follows. In (]2]we briefly recall the Boltzmann-type kinetic equation 
and we specialize it by giving two sets of interaction rules. The resulting S and x models are 
discussed in depth in ^3] and in © respectively. In particular, in ^3] we prove the existence of 
a class of quantized steady-state distributions. Since we are unable to prove their uniqueness, 
we discretize the model by approximating the kinetic distribution with a piecewise constant 
function and we then show by numerical evidence that the class of stationary solutions of the 
S model is only the one we have already studied analytically. In ^4] we show the somewhat 
surprising result that the equilibrium distributions of the \ model yield a macroscopic flow that 
is extremely well approximated by the discrete-velocity-based closure law resulting from the S 
model. This is illustrated in the final section £j5] where we show that the fundamental diagrams, 
that is the flux-density relationships, obtained from the two models tend to coincide under grid 
refinement. Next we compare these diagrams with experimental data, finding that our models 
reproduce well experimental fundamental diagrams and thus they capture the characteristics 
of macroscopic traffic flow. Further, we compute the macroscopic acceleration induced by the 
model, proving in particular its link with Av. We end the paper with a section summarizing the 
main results of this work, and proposing possible applications and further developments. Finally, 
the details of the matrix elements resulting from the discretization of the \ model are written in 
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an Appendix. 


2 The general form of the kinetic model 

In this section we briefly recall the general structure of a Boltzmann-type kinetic traffic model, 
which we will then specialize by prescribing a set of binary interaction rules in order to derive 
two models which differ only in the modeling of the acceleration interaction. Both models are 
defined on a continuous velocity space and they are characterized by a parameter Av related to 
the typical acceleration of a vehicle. 

We will focus on the space homogeneous case, because we want to investigate the structure 
of the collision term and of the resulting equilibrium distributions. In particular, we show that 
the simplified model ((5 model) permits to describe the complexity of the equilibrium solutions 
with a very small number of discrete velocities. 

Let / = f(t,v) : R + x V —> R + be the kinetic distribution function, where V = [0, Vmax] is 
the domain of the microscopic speeds and Vm ax is the maximum speed, which may depend on 
the mechanical characteristics of the vehicles, on imposed speed limits, environmental conditions 
(such as the quality of the road, the weather conditions, etc). The statistical distribution / is 
such that f{t , r))dw gives the number of vehicles with velocity in [u, v + dw] at time t. 

As usual, macroscopic quantities are obtained as moments of the distribution function / with 
respect to the velocity v: 

p(t ) = [ f{t,v)dv, ( pu)(t ) = [ vf(t,v)dv 

Jv Jv 

where p is the density, i.e. the number of vehicles per unit length (tipically, kilometers), u is 
the macroscopic speed and pu is the flux of vehicles. Note that p can also be interpreted as the 
reciprocal of the average distance between cars, see [2] ■ 

In the homogeneous case, the Boltzmann-type equation can be written as 

dtf(t,v) = Q[f,f](t,v) (1) 

where Q[f, f](t, v ) is the collisional operator which describes the relaxation to equilibrium due to 
the microscopic binary interactions among vehicles. For mass conservation to hold, the collision 
term must satisfy 

[ Q[fJ](t,v)dv = 0. 

Jv 

In fact, this ensures that, in the space homogeneous case, the density remains constant in time. 

The collisional operator is usually split into a gain term G[f, f] and a loss term L[f, /], that 
model statistically the interactions which lead to gain or to loose the test speed v. Denoting 
with A(z)*—>-u|u*) the probability that the velocity v £ V results from a microscopic interaction 
between candidate vehicles with velocity v* and field vehicles with speed v*, the model writes as 
an integro-differential equation 


dtf{t,v) = f [r](v*,v*)A(v^v\v*)f(t,v*)f(t,v*)dv*dv*-f(t,v) [ p{v,v*)f{t,v*)dv* (2) 
JvJv Jv 


G[fJ](t,v) 


L[f,f]{t,v) 


in which p(v*, v*) is the interaction rate possibly depending on the relative speed of the interacting 
vehicles, e.g. r/(u*,u*) = |u* — v*| as in HU [5]. Although such a choice would make the model 
richer, in |25| we found that a constant interaction rate is already sufficient to account for many 
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aspects of the complexity of traffic. Another possibility is to consider 77 as dependent on the local 
congestion of the road, that is r] = 77 ( 77 ). However this is not relevant in the homogeneous case, 
where p is constant, for then 77 would affect only the relaxation time towards equilibrium. Thus 
in this paper we will set 77 =constant. 

Notation. In the whole paper, in order to shorten formulas, we adopt the following traditional 
shorthand f(t, u*) = /*, f(t,v*) = /*, etc. Note in particular that in the space homogeneous 
case f* and /* are not different distribution functions, but the evaluation of the same f(t,v ) at 
two different points v* and v*. 

We will suppose that A depends also on the macroscopic density p in order to account for 
the influence of the macroscopic traffic conditions (local road congestion) on the microscopic 
interactions among vehicles, see (Ml EH HU ESI] ■ Thus, we suppose that A fulfills 
Assumption 1. 

A(u*— >v\v*] p) > 0, and / A(v*—tv\v*-,p)dv = l, for 77 *, v*, v € V, p G [0, p max ] 

Jv 

where p max is the maximum density of vehicles, for instance the maximum number of vehicles 
per unit length in bumper-to-bumper conditions. 

Remark 1. Any transition probability density A that satisfies Assumption [T] guarantees mass 
conservation since 

9 t [ f(t,v)dv= [ Q[f,f](t,v)dv = 

Jv Jv 

[ [ f(t,v*)f(t,v*)dv*dv*~ [ f(t,v)dv [ f(t,v*)dv* = 0. 

Jv Jv Jv Jv 

2.1 Choice of the probability density A 

The probability density A assignes a post-interaction speed in a non-deterministic way, consis¬ 
tently with the intrinsic stochasticity of the drivers’ behavior. The construction of A is at the 
core of a kinetic model. Here, it is obtained with a very small set of rules. 

• If 7 ;* < v*, i.e. the candidate vehicle is slower than the field vehicle, the post-interaction 
rules are: 

Do nothing: the candidate vehicle keeps its pre-interaction speed with probability 1 —Pi, 
thus v = i>*; 

Accelerate: the candidate vehicle accelerates to a velocity v > u* with probability Pi. 

• If 77 * > v*, i.e. the candidate vehicle is faster than the field vehicle, the post-interaction 
rules are: 

Accelerate: in order to overtake the leading vehicle, the candidate vehicle accelerates to 
a velocity v > 7;* with probability P2; 

Brake: the candidate vehicle decelerates to the velocity v = v* with probability 1 — P 2 , 
thus following the leading vehicle. 

From the previous rules, we observe that the probability density A has a term which will be 
proportional to a Dirac delta function at v = 77 *, due to the interaction which preserves the pre¬ 
interaction microscopic speed (the “Do nothing” alternative). Note that this is a “false gain” for 
the distribution /, because the number of vehicles with speed v is not altered by this interaction. 
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In the following, we assign the speed after braking as proposed in [22] and used also in BE] in 
the context of a discrete velocity model. Namely, we suppose that if a vehicle brakes, interacting 
with a slower vehicle, it slows down to the speed v* of the leading vehicle. Thus, after the 
interaction it gets the speed v = v* without overtaking the leading field vehicle. Instead, for the 
post-interaction speed due to acceleration we propose two different models. 

Quantized acceleration (<5 model): the output velocity v is obtained by accelerating instan¬ 
taneously from u* to the velocity min {w* + Av, U max }. Considering all possible outcomes, 
the resulting probability distribution, in this case, is 


A(v*->v\v*-, p) 


(1 T (u), if ft v 

(1 - P 2 )S V . ( V ) + P2S m in{ v ,+/\vy^}{v), if v* > v*. 


(3) 


Uniformly distributed acceleration (y model): the new velocity v is uniformly distributed 
between u* and min{i;* + Au, U max }. On the whole, the resulting probability distribution 
becomes 


A(v*—tv\v*] p) 


(1 — Pi)S v ,(v) + Pi 
(1-P 2 )6 v .(v) + P 2 


X[v* +Aji,V max }] ( v ) 

min{i)»+An,y m ax}—«» ’ 

X[xj»,miii{Tj„+ATj,V ma x}](' L ’) 

minlwj+AD.ynax} — V* ’ 


if < v* 
if f* > v*. 


(4) 


Note that the acceleration of a vehicle in ^ is similar to the one assumed in Big, which 
however were based on a discrete velocity space. In HI IS] the acceleration parameter Av is 
chosen as the distance between two adjacent discrete velocities, thus Av depends on the number 
of elements in the speed lattice. In this work, Av is a physical parameter that represents the 
ability of a vehicle to change its pre-interaction speed w*. With this choice, Av does not depend 
on the discretization of the velocity space and the maximum acceleration is bounded, as in HU. 
In contrast, deceleration can be larger than Av, and this fact reflects the hypothesis that drivers 
tend to brake immediately if the flow becomes more congested, while they react more slowly 
when they can accelerate (see the concept of traffic hysteresis in |28j and references therein). 

The acceleration performed in the y model has some points of contact with the microscopic 
rules prescribed in [T2i 114] . In [12], however, the post-interaction speed is selected through a 
random process in the interval \v* , V max ]. Here instead, the post-interaction speed is determinis¬ 
tic. In m instead, the velocity after acceleration is uniformly distributed over a range of speeds 
between u* and v* + a(V me , x — v*), where a is supposed to depend on the local density; in a similar 
way, the output velocity from a braking interaction is assumed to be uniformly distributed in 
[/3v*,v*], with /3 £ [0,1]. 

In the following, the probabilities P± and P 2 are taken as P\ = P 2 =: P and P will be a 
function of the local density only, as assumed for instance in [23] where P = 1 — p/p max . More 
generally, from a modeling point of view, P should be a decreasing function of p, see also m or 
|2!|- For instance in [23] we have considered 


P= 1- 


(5) 


where 7 £ (0,1) can be chosen to better fit experimental data. In [7] the choice P = a(l—p/p max ) 
is proposed, a £ [ 0 , 1 ] is a parameter describing environmental conditions, for instance road or 
weather conditions. In more sophisticated models, one could also choose P as a function of the 
relative speed of interacting vehicles, but we will not explore this possibility in the present work. 

The simplified choice Pi = P 2 and the interaction rules described at the beginning of this 
section guarantee the continuity of the transition probability ([3| and Q along u* = v*. 
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Remark 2. Both choices © and © for A include terms of the form 5 Vt (v), which actually 
describe false gains mentioned above, because the velocity of the candidate vehicle does not 
change. They are automatically compensated by false losses, as it can be seen by rewriting the 
classical kinetic loss term of equation © in the form 

L [f,f](t,v)= [ [ r]S v „(v)f*f*dv*dv*. 

Jv Jv 


3 The (5 velocity model 


Now, we focus on the steady states of model ([3]). We start with the existence of a particular 
set of equilibrium solutions of the continuous model, which are computed analytically. Next, we 
consider a finite volume discretization of the model, and we show that the discrete equilibria 
have precisely the structure found before analytically, thus suggesting that the particular set of 
equilibria found analytically are the only equilibria of the system. 

It can be proven |S] that the Cauchy problem associated to ([2]) is well posed provided the 
probability density A is Lipschitz continuous with respect to v* and v* in a suitable Wasserstein 
metric. This is indeed the case of the A defined in © , with Pi = P 2 . 

Using the expression © for A, we rewrite the gain term in © as 


=V [ [ [(! - P)6mm{v.,V}(v) + f , £ m in{t,,,+Au,V max } («)] f*f*dv*dv* 

Jv Jv 

and the following important result on the existence of a particular class of stationary solutions 
holds. More precisely, it characterizes the equilibrium distributions having the form of linear 
combinations of Dirac’s masses. This theorem establishes the connection represented by the left 
vertical arrow in Figure [T] 


Theorem 3. Let P be a given function of the density p G [0,p max ] such that P G [0,1]. Let 
{ v j}jLi be a set of velocities in [0, Umax]- The distribution function 


N 

/”(») = Et 

i=i 


°°A 

o 


». /r>° v i = i > 


,JV, 


with l fj° = Pi is a weak stationary solution of the S model provided Vj = V\ A jAv, 


Proof. Without loss of generality suppose that is an ordered set of velocities such that 

0 < V\ < ■ ■ ■ < Vn < U max . Since the distribution function f°°(v) = Y^f=i fj°^ v A v ) a wea k 
stationary solution of the S model, it satisfies the following steady weak form of equation ©: 


f f ( f (f^Afv^-tvlv*-, p)dv\ f 00 (v*)f°°(v*)dv*dv* - p [ <p(v)f°°(v) du = 0, 

Jv Jv \Jv J Jv 


where (j) € C c {V) is a test function, with C c (V) the space of continuous functions having compact 
support contained in V, and the probability density A is given in ©• Substituting the expression 
of f°° in the above equation we obtain 


N k N N 

f 1 - p ) + ( 1 -- p )£ E 

k= 1 h—1 k—1 h=k -\-1 

N N 

+-Pp E <K min {^+ Av , Kna x})/r - p E msr = °- 

h= 1 3 =1 


( 6 ) 
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The proof will be organized as follows: in order to determine an equation for the /?°’s, we 
consider a particular family of test functions (p 7 defined as piecewise linear functions such that 
4>j{vj) = 1 and 4>j(vi) = 0, V i A j; in this way, first we find an equation for /f°, then we show 
that /£° 0 if v -2 = Ui + Au and finally by induction we prove that if f°° A 0, then Vj = Vi+jAv, 

for some j £ {1,..., N}. 

Let j = 1, equation Q with <j> = (f >i becomes 

N N N 

a /r/r + (! - p) E fhf? 0 +ppJ2^ fa+ Av )/“ - pf?° = °- 

fc= 1 h =2 h=l 

Due to the particular construction of fa and using Yif=i fj° = Pj the above expression reduces 
to 

-(i-p)(/n 2 + (i-2P)p/r = o 

which admits the two roots //° = 0 and /f 3 = p 1 1 J 2 / > • If P > 1/2, only /^° = 0 is acceptable, 
because the other root is negative. If instead P < 1/2, both roots can be accepted, but only 
/f° = Pqpp- > 0 is stable. This argument will be used for selecting a single root throughout the 
proof. 

Now, let j = 2 and <j> = <j > 2 . Equation © writes as 

N 

- (1 - W 2 00 ) 2 + [(1 ^ 2P)p - 2(1 - P)A°°] f? + Pp^ ^(min{« fc + Av, V ma xM°° = 0. (7) 

h =1 


Note that fa (min{ Vh + Au,I/ max }) = 1 if h = 1 and Vi + Av = V 2 ■ While ^(min \yh + 
Av , A na x}) = 0 otherwise, due to the particular choice of fa which is centered in w 2 and can be 
taken with support smaller than 2 |w 2 — v\ — Aw] around this point. Then, if v 2 A V\ + Av , the 
constant coefficient of ([7]) is zero for all P £ [0,1] and the solutions of the equation are f%° = 0 
or / 2 0 = p 2 p — 2/f°. Exploiting the structure of /f° and the fact that /°° > 0, only /(p = 0 
is the admissible solution for all values of P £ [0,1]. 

Instead, if u 2 = tq + At:, the third term in equation ([T]) is Pp/f° > 0, for all P £ [0,1]. 
More precisely, if P > 1/2 then /f° = 0, thus the constant coefficient vanishes and again 
one concludes that /|° = 0. While if P < 1/2, Ppf£° is positive and since the discriminant 
V = ((1 — 2 P)p — 2(1 — P)//°) 2 + 4P(1 — P)p/i° of equation ([T]) is positive and the leading 
coefficient is negative, the equation has two real roots with opposite signs. Therefore 


/ 2 ° 


-(1 - 2 P)p + 2(1 - P)f f° - fav 
—2(1 — P) 


and this is the only case in which f%° can be non-zero. 

We now proceed by induction. Suppose that Vk — v± is an integer multiple of Aw and 


fk 


°, 

-(l-2P)p+2(l-P) fr-VTh; 

— 2(1 —P) 


if P > 1/2 
if P < 1/2 


for all k = 3,... ,j - 1, where V k = ((1 - 2 P)p - 2(1 - P) Yti fi°) ' + 4P ( X - P )p YiZl fi° ■ 
We show that f°° can be non-zero only if Vj = tq + j Av, for j £ {1,..., N}. 

Taking the test function (j) = the equation for f°° writes as 


N 

fj° + PpYM^Vh + Al ’> Kiax})/r = 0. 

h =1 

( 8 ) 
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(1 -p)UT) + 


i -1 


(l-2P)p-2(l-P)^/r 


Z=1 









Note that <j)j + Av, Kn ax }) = 1 if h = j — 1 and Vj-± + Av = Vj. While (f>j(mm{vh + 
Av, Vmax}) = 0 otherwise, due to the particular choice of 4>j which is centered in Vj and can be 
taken with support smaller than 2 | Vj — Vj —i — Ai>| around this point. If Vj = V\ + jAv, then 
Vj = Vj-i + Av and the above equation becomes 


-(1 -p){frr + 


j -1 


(l-2P)p-2(l-P)^/f 


1=1 


fj° + PpfV -1 = 0 


which has two real roots. If P > 1/2, using the inductive hypothesis /£° = 0 for all k < j — 1. 
Thus f°° = p^jzrpr (which is again not acceptable since it is negative) or f°° = 0, confirming the 
induction. If P < 1/2, we have two real roots with opposite signs, so one proves that f°° can be 
chosen strictly positive. 

If instead v 3 ^ v\ + j Av then the constant term of equation ([8]) is zero and the two roots are 

fj° = 0 and f°° = Sj = l 1 2P ' IP tl which is negative for all values of P £ [0,1] 

using the Lemma [4] below. □ 

In the previous proof we use the following technical fact. 

Lemma 4. Let P be a given function of the density p £ [0, p max ] such that P £ [0,1]. Consider 
{.fj}jLi € defined as 

, _ 2(1-P)ECi/^(1-2P)p-C j 

h -2(1 - P) 

with Cj > 0. Assume that fj is positive for all j. Then 

(l-2P)p-2(l-P)Et- 1 1 / i _ 

= - —- - < 0 


for all 2 < k < K. 

Proof. We proceed by induction on k. Let k = 2 

c (1 — 2P)p — 2(1 — P)fi C\ 

52 = -T^p- = ^i^p < °' 

Suppose Sj < 0, j = 3,. .., k — 1, then 5/ < 0. In fact 

S k = S k . i - 2/*_! <0. □ 


3.1 Discretization of the model 

Observe that Theorem [3] ensures the existence of a class of steady solutions for the <5 model 
which are characterized by the fact that the total mass of vehicles on the road is distributed only 
on the velocities which are multiples of Av. We cannot prove the uniqueness of such a class of 
steady solutions. However we can show by numerical evidence that the asymptotic distributions 
of the S model are only of the type stated by Theorem [3j Thus, in this subsection we introduce 
a discretization of the model ([3|. 

To this end, the explicit formulation of the gain term is now useful. Notice that the Dirac 
delta function at v = min {u* + Av , V mSiX } can be split as 




(^)i if V* £ [O5 f/nax Au] 

if ^ (Imax ArJ, 1/nax] 
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because the velocity jump of size Av, leading to the output velocity v = v*+Av, can be performed 
only if u* < Vmax ^ Av. If instead u* £ (Vmax — Av, Vmax], the post-interaction velocity will be 
v = Vmax. Thus the gain term of the S model can be written as 


G[f, /](*, v) = 77(1 - P)f(t , v ) 




/*du* 


+ rjPp 


f(t,v 


A v)H Av (v) + 6 Vm „{v) 



(9) 


where H a {x) denotes the Heaviside step function with jump located in a. The last term in the ex¬ 
pression of G means that, as a result of the microscopic interactions, the mass Pp f ” ax * /*dv* 
is allocated entirely to the speed Vmax- Note that, in space-nonhomogeneous models, /» and /* 
may refer to distributions evaluated at different locations in space, see for instance m and m- 
For this reason we keep the integrals over field and candidate particles separate. 

Suppose for simplicity that the acceleration parameter Av satisfies Av = V max /T with T £ N. 
We consider a discretization of the velocity space defining the velocity cells Ij = [(j — (j — 
D [0, Vmax], for j = 1,..., N. Note that all cells have amplitude Sv = Vmax/ {N — 1) except 
1 1 and In which have amplitude 5v/2. 

We consider a piecewise constant approximation of the kinetic distribution so that 

N ( x 

f(t,v )« f N {t,v) = (10) 

i=1 

where fj represents the number of vehicles traveling with velocity v £ Ij. 

By integrating the kinetic equation ([I]) over the cells Ij and using /jv(t,v) in place of f(t,v ) 
we obtain the following system of ordinary differential equations 


fj(t) = Qj[f,f](t) ■■= [ Q[f,f](t,v)dv, 
Ji 0 


( 11 ) 


whose initial conditions /i(0),..., /at(0) are such that 



f(t = 0,u)du = p 


and p is the initial density, which remains constant during the time evolution in the spatially 
homogeneous case. 

We set r := Av/5v £ R + and we define r + := r+ 1, r~ := r — 1. Then /[>+] is the cell which 
contains v = Av , where |"r + ] denotes the integer part of r + . By computing the right hand side 
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A 3 S for 1 < j < \r ] A 3 S for j = |Y + ] 



Figure 2: Structure of the probability matrices of the S model, with Av = Vmax/2 













of the ODE system (11), we obtain explicitly 


N 


1 k=j +1 


(12a) 


N 


N 


a -p)f j E A-/iE/ fe ’ 


for i = 1,..., |r ] 


h=? + l 


fc=l 


TV 


N 


-Q j {fj](t) = (i-p)f j f j +(i-p)f j E / fc + (i-nr E a 


(12b) 


fc=j+l 


h=J+l 


+ Pp 

-/iE/ 


2/imin{i, [r+] - ^ - r} + ,fr—T/ 2 (E 1 “ r ) 


AT 


for j = |~r + ] 


fc=i 


AT 


AT 


-Qi[/./](*) = (i-W/i + (i--P)/i E / fc + (i-W E A 

^ fe=i+i ft=i+i 

+ Pp [(1 + < ^H,[r+l^,|'r+l+l)/i-|'rl (1 + r EEltjEjitEZljt [ifl! 

N 

-/jE^’ for j= |r+] +l,...,iV-l 


(12c) 


k =1 


: Qjv[/,/]w = (i-p)r/ J v 


(12d) 


+ -Pp^M.rr+l 

-Pp/iV-Tr-l 
N 

-f N J2f k - 


/at—|Y+] (r — \r l) + /iv-rr-l([r + ] - - - r) 


+ 0- V 1 + \) 


N 


p P E A 

h=N— [>+]+2 


fc=l 


where here <5 t j’s are Kronecker’s delta’s. The terms with a wavy underline are those deriving 
from the acceleration term. In the formulae above, the position of the index of the components 
of f = /at] t € distinguishes the distribution of the field and of the candidate vehicles: 

bottom right for the candidate vehicles (as in fh ), top right for the field vehicles (as in f k ). In 
vector form: 


_d 

dt 


fj = V 


f T ^f - f T e,lT f 




(13) 


1* 


where denotes the vector with a 1 in the jr-th component and 0’s elsewhere, 

[1,..., 1] £ I w . The matrices A J s have a sparse structure, shown in Fig. [ 2 ] in which the nonzero 
elements are shaded with different hatchings, corresponding to the different values of the ele¬ 
ments, as indicated in the panels in which they appear for the first time. 

As it can be checked using (121, these matrices are stochastic with respect to the index j, 
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Parameter 

Description 

Definition 

N 

Sv 

r 

T 

number of discrete speeds 
cell amplitude 

ratio between the speed jump Av and the cell amplitude Sv 
number of speed jumps Av contained in [0, V ma , x ] 

Sv = 

r 

1 

N-l 

Av 

y v 

Av 


Table 1: Table of the numerical parameters. 





Figure 3: Structure of the probability matrices of the 6 model with Sv integer sub-multiple of 
Av. 


i.e. £f=i = 1, V/i, k e {1 This property comes from Assumption jlj and it 

guarantees mass conservation. 

Recall that the elements of the matrix are the probabilities that the candidate vehicle 

with velocity in Ih, interacting with a field vehicle with velocity in Ik, acquires a velocity in Ij. 
The fact that these matrices are sparse means that a velocity in Ij can be acquired only for 
special values of the velocity of candidate and field vehicles. In particular, the j- th row of the 
matrix A^ contains the probability of what we called “false gains” in Remark ij that is the 
probability that the candidate vehicle does not change its speed. The non zero elements of the 
j- th column are the probabilities that a candidate vehicle acquires a speed in Ij by braking down 
to the speed of the leading vehicle. The non zero rows, located at h = j — [r"| |>+], h — 1 and 
h + 1, contain the probabilities that the candidate vehicle accelerates by Av, acquiring therefore 
a velocity in u* + Av £ Ij starting from a velocity u* in Ih-i, Ih or Ih+i- The band between the 
rows j — |V| + . r r +i + 2 and j — 1 is filled with zeros, because in the <5 model the acceleration 

is quantized. As we will see in Section [4.1[ this band will be filled by non zero elements in the \ 
model, where the acceleration is distributed uniformly between [0, An]. 

In Table [I] we summarize the numerical parameters introduced in order to discretize the 
continuous-velocity model. 


Sv as integer sub-multiple of Av. For the special choice of the velocity grid which ensures 
that r = Av/Sv £ N, then the formulae (121 simplify and the resulting interaction matrices are 
given in Figure [3] Notice that the rows j — r ± 1 are filled with zeros. In fact, since Sv is an 
integer sub-multiple of An, a velocity in Ij can be obtained as result of an acceleration only if 
the pre-interaction speed is a velocity in Ij— r - 
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Figure 4: Approximation of the asymptotic kinetic distribution function obtained with two 
acceleration terms Aw = 1 /T, T = 3 (top), T = 5 (bottom), and N = rT + 1 velocity cells, with 
r € {1,4,8}; p = 0.3 (left) and p = 0.6 (right) are the initial densities. We mark with red circles 
on the x-axes the center of the T + 1 cells obtained with r = 1. 


The structure of the matrices A J S determines the equilibrium of the discrete model (131. In 
Figure |4]we show the function (w) = Hindoo /)v(t,w) obtained by integrating numerically the 
system of equations up to steady state, for a few typical cases. In all numerical tests we take 
P = 1 — p, with V max = p m ax = 1- As initial macroscopic densities, we choose p = 0.3, 0.6 
(plots to the left and right of the figure). We consider two values for the acceleration parameter, 
Av = Fmax/T, T = 3, 5 (top and bottom of the figure). The number of velocities in the grid is 
taken as N = rT + 1, with r £ {1,4, 8}. The three curves in each plot contain the data for the cell 
averages of the equilibrium distribution for the different values of r. It is clear that in all cases, 
/jy (w) is a function of the density p and it approaches a series of delta functions, centered in the 
velocities which are multiples of Aw and indicated in the picture by red dots on the horizontal 
axis. 

This means that, as 5v — > 0, only a finite number T + 1 of velocities carry a non-zero mass of 
vehicles at equilibrium. More precisely, the discrete asymptotic function /“(w) is different from 
zero only in the T + 1 cells I\, I r , / 2 r ,..., In- Therefore, as time goes to infinity the number of 
nonzero values of the fj's appearing in (10) is univocally determined by the acceleration term 
Aw = V ma x/T. 

The previous considerations can be supported by the numerical results in Figure [5j in which 
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Figure 5: Evolution towards equilibrium of the discretized model (131 with TV = 4 (green), TV = 7 
(blue) and TV = 10 (red) grid points. The acceleration parameter Av is taken as V max /3 and the 
density is p = 0.6. Black circles indicate the equilibrium values. 


we show the evolution towards equilibrium of the fj' s, j = 1,..., TV. In this figure, Av = V max /3, 
and the different plots are obtained starting from a uniform initial distribution, namely fj(t = 
0) = p/N,\/j = 1,..., TV, with r = 1 (green), r = 2 (blue) and r = 3 (red), which correspond to 
TV = 4, 7 and 10 velocity cells respectively, and p = 0.6. It is clear that under grid refinement 
the number of nonzero steady values does not change. In fact, note that a different dynamics 
towards equilibrium is observed, for different values of the number TV of cells, but as equilibrium 
is approached, the values of the fj' s go to zero except for the velocities corresponding to integer 
multiples of Av. Moreover, the non zero values of the steady-state distribution /°° do not depend 


on the discretization parameter Sv. This fact can be also deduced by looking at equations (12) of 


the discrete collision operator. These expressions are not functions of Sv. Thus all exact values 
of the equilibria can be obtained using a coarse grid. 

Theorem [6] below confirms the structure of the equilibria that we have just observed in the 
numerical results and it states that the steady-state solution of the S model, prescribed by 
Theorem [3j can be reconstructed numerically on the grid with Sv = Av. To this end, we first 
recall a result from |8J, where the existence and well posedness of the solution of such systems is 
proved and then we show that the discrete equilibria of additional equations resulting from the 
choice r > 1 give actually no contribution. 

Theorem 5. (Delitala-Tosin) Let fj(t = 0) > 0, with )Cyli fj(t = 0) = p, be the initial condition 
for the system 

d -Jj = t T A>f - f T ej l T N f, j = 1,.... TV, 


dt“ 


where the matrices are stochastic matrices with respect to the index j, i.e. Yhj A J hk = 1 for 
all h, k. Then there exists f* = +00 such that the system admits a unique non-negative local 
solution / > 0 satisfying the a priori estimate 

\\f(t)\\i = II/0II1 = P VT G ( 0 ,t* = +00]. 

The following result, together with Remark [TJ shows that all equilibria of the discrete model 
are of the quantized form described by Theorem [3j Thus the next Theorem establishes the 
correspondences symbolized by the right vertical and the middle horizontal arrows in Figure |T| 
Theorem 6. Let P be a given function of the density p. For any fixed Av = V max /T, T £ N, let 
f r (p) denote the equilibrium distribution function of the ODE system (131, obtained on the grid 
with spacing Sv given by Av = rSv with r = (TV — 1)/T G N. Then 


(fr),- = 


(fi) |"ij mod(j — 1, r) = 0 


0 


otherwise 
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is the unique stable equilibrium and the values of fi depend uniquely on the initial density p , 

with Efc=i( f i)fc = P- 

Proof. We already know from Theorem [5] that the solution of (131 exists, is non-negative and is 
uniquely determined by the initial condition. 

To prove the statement, we compute explicitly the equilibrium solutions of (131, using the 
explicit expression of the collision kernel given in (12aI, (12b), (12c) and (12d), with r £ N and 
N = rT + 1. Since here we are interested in the solutions of the homogeneous problem, we will 
take identical distributions for the candidate and the field vehicles, i.e. fj = f'P 


For j = 1, using the expression (12a I and the fact that Efc=i fk = P we obtain 

-(1 - P)fl + (1 - 2P) pfi = 0. 




<£> 


This is a quadratic equation for /,, which has the two roots /i = 0 and /, = p( 1 — 2P)/(1 — P). 
It is easy to see that one solution is stable, and the other one unstable, depending on the value 
of P. Here we are interested only in the stable root, so we find 


(fr)l = 


0 

Ppf 


P > - 
r — 2 

otherwise 


(14) 


Thus, no vehicle is in the lowest speed cell I\ if P > which, for the simple case P = 1 — p 
means that all cars are moving if p < 

The case j = 1 we just computed is typical. Also for larger values of j, we find a quadratic 
equation for the unknown fj, which involves only previously computed values of fk , k < j. Thus, 
we can easily compute iteratively all components of f r . 

For 2 < j < r, the equilibrium equation is obtained by using again the expressions (12a I with 
r £ N. Then 

I" j—l 

-(1 -P)/? + (1 - 2P)p - 2(1 - P) £ f k fj = 0. 

k =1 . 

Start from j = 2. Clearly, for P > substituting equation (14), we again have (f r )2 = 0. For 
P < \ , the equation for f 2 , with /, given by (141, becomes 

-(1-P)/|-(1-2P) pf 2 = 0. 

Comparing with the equation for /i, we see that now the stable root is f 2 = 0. Thus, at 
equilibrium, we have (f r ) 2 = 0, for all values of P. Analogously, it is easy to see that (f r )j = 0, 
Vj = 3,...,r. 

For r +1 < j < N, in place of (12a), we use (12b) and (12c I in order to obtain the equilibrium 
equation. Then 


- (1 - P)fl + 


3 ~ 1 


(1 - 2 P)p - 2(1 -P)J2 fk 


k = 1 


fj + Ppfj-r = 0. 


The equation has a positive discriminant 


T>i = 


3 -1 


(l-2P)p-2(l-P)^/ fc 


k = 1 J 


+ 4P(1 — P)pfj_ r 


( 15 ) 
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thus it always admits two real roots. To fix ideas, let us consider r + 1 < j < 2r. If j = r + 1, 
since (f r )k = 0, Vk = 2,,r, equation (151 becomes 

-(1 - P)f 2 r+1 + [(1 - 2 P)p - 2(1 - P)h] f r+1 + Ppf x = 0. 


Thus, we find (f r ) r +i = 0 for P > 1/2, because the equation for f r+1 becomes identical to (3.11. 
If instead P < 1/2, substituting the expression for /i, the equation for f r+ 1 becomes 


-(l-P)/ r 2 +1 -(l-2P)p/ r+1 + p 


2 P(1-2P) 


= 0. 


r+ l , r 1 _ p 

This equation has a negative and a positive real root, which is stable. Thus 


[° 

(f r ) r+ l — ^ _(i_2P) p + Pa /(1-4P 2 ) 

l 2(1 —P) 


P >-2 

otherwise. 


Now, let r + 1 < j < 2r. Since fj- r = 0, the constant term of equation (151 is zero. Then, as 
seen for j = 2,... , r, it is easy to prove that, for each j = r + 2,... ,2r and for all values of P, 
one solution is negative and thus (f r )j = 0. 

Clearly, this procedure can be repeated and we find that 


0 


(fr)i - < _ 2 (l-P)]T J fc= 1 i / fc + (l-2P)p+7p“ 


2(1—P) 


P > i 
r — 2 

otherwise 


if j = Ir + 1, l = 0,..., T — 1, while (f r )j = 0 otherwise. From these considerations, the thesis 
easily follows. Finally, just note that for the last value, /n, we can use mass conservation 


T -1 

(fr)rT+l = (fr)lV = P — ^ (fr)ir+l (p) ■ 
1=0 


□ 


Notice that the result above proves that equilibria are determined by the initial density p, 
which is constant in time in the spatially homogeneous case, but they do not depend on the 
number of cells N used to approximate the kinetic distribution. In fact, the stable equilibria just 
computed correspond exactly to the values f°° given in Theorem |3] and they can all be recovered 
on the coarse grid Sv = Av, i.e. choosing r = 1. 


Remark 7 (The case of a generic Sv). In order to further investigate the existence of stable 
equilibria, we can seek more general ones numerically. In fact, the finite volume discretization 
(10 ) is capable of converging to absolutely continuous equilibria, but so far we proved that it 
converges to sums of Dirac masses if Sv is an integer submultiple of Av. Here we apply our 
discretization scheme with non-integer ratios r = Av/5v and show that the resulting equilibria 
converge, in the sense of distributions, to the equilibria already described in Theorems [3] and |6] 
When r = Av/5v is not integer, Theorem [h] cannot be applied, but numerical integration 
of equation (131 shows that, for large time, the solution approaches equilibria that have masses 
concentrated at points spaced by Av. More precisely, in this more general case, only a small 
finite number of f °°'s is nonzero and the cumulative distribution function induced by the discrete 
/jv», cf- pk, 


pV 

F n (v)= / /n(v) du, J)£[0,k m 
Jo 
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p=0.6, T=3 


p=0.6, T=5 




Figure 6: Cumulative density at equilibrium for several values of Sv —> 0. The density is p = 0.6 
and Av is chosen as 1/3 (left), 1/5 (right). 


approximates the cumulative distribution of a sum of Dirac masses centered at multiples of Av. 
That is, Fpf(v) converges to a piecewise constant function with jump discontinuities at multiples 
of Av. See Figure [bj 

Furthermore, Figure [6] shows that the jump discontinuities of the cumulative distribution 
in the limit Sv —> 0 are located exactly in the points computed analytically by Theorem [6] 
in the case r £ N. In particular, the figure shows the cumulative distribution at equilibrium 
computed by solving numerically (131 for several values of the discretization parameter Sv with 
non-integer ratio A v/Sv. In both panels the density is p = 0.6, while Av = 1/3 in the left plot 
and Av = 1/5 in the right one. The continuous red curve represents the cumulative distribution 
of the stationary solution for Sv = Av, computed by Theorem [6j The black dotted curve is 
computed with N = 15 velocity cells, the green clot-dashed line with N = 30 and finally the blue 
dashed one with N = 60. We observe that as Sv —> 0 the cumulative distribution tends to that 
of the analytic solution of Theorem [3] where only the velocities centered in multiples of Av give 
rise to a jump discontinuity. 


We conclude the section with a few remarks on the structure of the equilibria of (13 ). 


Remark 8 (Reduced velocity space). The S model is characterized by a small number of non¬ 
trivial values for the microscopic velocities, which allows one to compute analytically the equi¬ 
librium of the system. This fact provides an analytic closure of the macroscopic equations. This 
fact can also be exploited from a numerical point of view, justifying the use of a small number 
of discretization points in simulations aimed at capturing equilibrium effects. In this sense, the 
model seems to suggest that kinetic corrections can be accounted for with a computational cost 
which is not much higher than the one needed for a macroscopic model. 


Remark 9 (Unstable equilibria). Theorem [6] gives the uniqueness of the stable equilibrium of 
the model with Av/Sv = r £ N. Unstable ones may occur if the initial condition is such that 
/i(0) = 0. In fact, the interaction rules related to the case v* > v* do not generate a post¬ 
interaction velocity v which is less than v*. Thus if /i(0) = 0, i.e. if there are no vehicles with 
velocity v\ = 0 at the initial time, there will not be interactions leading to an increase of f±. This 
consideration can be generalized: if fj( 0) = 0 for j = 1, ...,/< r, then the computed equilibria 
will be fj = (f r)j_ji where f r is the vector containing the stable equilibria. In this sense, the 
equilibrium solution of the S model does not only depend on p but also on the initial condition 
f(0,v). These solutions however are unstable: a small perturbation on f\ (t = 0) is enough to 
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Figure 7: Evolution towards equilibrium, p = 0.7, T = 4, N = 17. Left: fj(t = 0) = p/N. 
Middle: f :l (t = 0) = 0 ,j = 1,2,3, J)(t = 0) = (p/(IV - 3)), j > 3. Right: /, = e = lO” 6 , / 2 = 
/3 = 0 and fj(t = 0) = ((p — e)/(iV — 3)). The thick lines highlight the components fj and the 
blue ones are for those that appear in stable equilibria, i.e. with j = kr + 1 for k = 0,... ,T. 




Figure 8: Speed of convergence towards the stable equilibria of the 8 model. The initial condition 
is a small random perturbation of the steady-states. 


trigger the evolution towards the stable equilibrium, which depends only on p. 

This is illustrated in Figure [7] in the left panel we show the evolution towards equilibrium 
when fj (t = 0) 0 for all classes (this is the stable equilibrium), while in the middle we show the 

case when fi = f 2 = f 3 = 0. In the rightmost panel we show a perturbation of the previous case, 
where /i takes a very small but nonzero value. It is clear that the evolution goes at first towards 
the unstable equilibrium of the middle panel, but then, in the long run, the stable equilibrium 
of Theorem [G] emerges. 


Remark 10 (Convergence rate to equilibrium). In Figure [8] we study the rate of convergence to¬ 
wards the stable steady-states. For the set of densities p £ {0.2, 0.3,0.4, 0.6,0.7, 0.8} we integrate 
numerically the system (13) for large times starting from a small random perturbation of the sta¬ 
ble equilibrium. In both panels, we use a linear scale for the x-axis (time) and a logarithmic scale 
for the y-axis (error). The error is computed at each step as e(t) = || f{t) — /°°|| 2 . The figure 
suggests that the rate of convergence towards the stable equilibrium depends on the density. In 
fact, in the left panel we consider different values of the ratio A v/8v = r, in particular r = 1 (solid 
lines) and r = 2 (dashed lines), and we note that the slopes of the curves corresponding to the 
same value of p are the same. A similar behavior is observed in the right plot in which two differ¬ 
ent Av are taken, Aw = 1/5 (solid lines) and Aw = 1/3 (dashed lines). We can conjecture that for 
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large enough times the distance from equilibrium behaves as e{t) ~ Ce M G) t \\f(t = 0) — /°°|| 2 , 
where C = C(r,Av), M(p ) > 0. 


4 The x velocity model 


The structure of the steady-state distribution of the S model clearly depends on the particular 
choice of the acceleration interaction made in ([3]), in which a vehicle accelerates by jumping from 
its pre-interaction velocity u* to the new velocity v* + Av. Thus it could seem quite natural that 
only velocities 0, Av, 2Av,..., V max give a non zero contribution at equilibrium. 

Here we study the x model, already introduced in Section 2.1 (see equation Q), in which 
vehicles can assume a post-interaction velocity uniformly distributed over a range of speeds when 
the acceleration interaction occurs. We will show that, although this model is more refined than 
the S model, at equilibrium the essential information is already caught by the simpler S model. 

Using the formulation Q for the transition probability density A, we rewrite the gain term 
in ([ 2 ]) as 


— V* 


G[fJ](t,v) = r 1 f [ (1 + 

JvJvl mm {v* + Av, V max } 

Notice that the % function can be split as 

X[t)«,min{u,-|-Ajj,y max }] { v ) j ~ ~ ~ > if £ [0, Vjnax — Ay] 


f*f*dv*dv*. 


min (u* + Av, U max } - t>* 


X[v 


1 ( V) 
XJ v > 


V max -u* 


if V* G (^max A'U, V^nax] 


hence substituting in the above equation and evaluating explicitly the integrals, we find 

C Knax r Knax 


G[f, f] (t, v) = 77(1 - P)f{t, v) 
+ r)Pp ' 


f*dv* 


f*dv* 


Vmax-At; 


X[v„v.+Av](v)f*dv* 


rV m 
Vnax — 


X[t;»,V max ](' 1 ’) 


A I /VI 'O'* -f-AAC' 1 I x r 

Jo J Vmax-At; Vmax ~ V* 

Observe that (16) differs from © only in the terms proportional to P. 


f* du» 


(16) 


4.1 Discretization of the model 

To compute the steady-state solution of the x model, we need to integrate the equations nu¬ 
merically. We use the same discretization of the velocity space V introduced to discretize the S 


model and therefore the kinetic distribution is approximated as in ( 101 . 

Integrating the kinetic equation (Jl|) over each cell, we find the system of ODEs (111, but 
now the gain term is given by ©• Although in this case the integrals are laborious, they 
can be computed recalling that Av = V max /T with T £ N and assuming that Sv is an integer 
submultiple of Av. Thus we will take N — 1 = 0 (modT) and r = 

The details are reported in Appendix [A] Here we just point out that, as in the case of the <5 
model, the ODE system can be conveniently rewritten in vector form as 

tJj = V - f T e,lT f] , j = 1,... ,N 


dt 


(17) 

where f = [f±, ..., /at] t £ R w is the vector of the unknown functions, ej £ denotes the 
vector with a 1 in the j-th coordinate and 0 ’s elsewhere, 1 ^ = [ 1 ,..., 1 ] £ R^ and A 3 X is the 
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N 

I 


3 ~r 


N — r 


N 


A 3 for 1 < j < r 


A\ for r < j < N 


A{ for j — N 


Figure 9: Structure of the probability matrices of the \ model with Sv integer sub-multiple of 

Av. 


j -th interaction matrix such that (A-Q contains the probabilities that a candidate vehicle with 
velocity in If, interacting with a field vehicle with velocity in I ^ acquires a velocity in Ij. 

Since the matrices appearing in are again stochastic, we can apply Theorem [5] to guar¬ 
antee the well posedness of the associated Cauchy problem. 

In Figure [9] we show the structure of the \ matrices: they are less sparse than the S matrices 
because of the uniformly distributed acceleration in [u*, u* + Au], where u* is the pre-interaction 
speed. In fact the matrix A 3 contains non-zero elements also in the rows from the (j — r + l)-th 
to the (j — l)-th, see the shaded areas in Figure [9] which represent non-zero probabilities of 
accelerating to a speed in Ij. Since Av = rSv, exactly r rows fill up. Instead, the area drawn 
using hatching contains the same probabilities already shown in Figure [3] for the case of the 
S model, with r G N. Note that the elements of the matrices depend on Sv (see equations 
(A.l l). Thus, in contrast to the <5 model, steady solutions of the ODE system ( [T7| depend on the 
number of velocity cells N chosen to approximate the kinetic distribution (see (fTo|). In other 


words, although Sv is an integer sub-multiple of Av, this model does not converge, as time goes 
to infinity, to the asymptotic distribution on the coarse grid as in the S model. However, we do 
recover the asymptotic distribution on the coarse grid in the limit Sv —> 0. 


A 3 
x' 

A j 

x 


Notice from equations (A.l I that all the elements in the rows j — r,... ,j — 1 of the matrices 
, j = 1 tend to 0 as 1/r when the grid is refined. In particular, for j = l,...,r, 

Ag. This consideration is not true for the matrices A 3 for j = r + 1,..., N. 


4.2 Expected speed of the 8 and the \ model 

Despite their differences the \ and the S model are deeply related. This can be seen by computing 
the expected output speed in each model resulting from a fixed pre-interaction speed. We define 
the expected value (v) of the post-interaction velocity as 


fV m 


(v) = 


vA(v*—>-vjv*; p) dv, (v*,v*) G V x V. 


(18) 


For brevity we indicate with A$(v) and A x {v) the probability densities given in (|3j) and Q 
respectively. Again we assume that Pi = P 2 = P as function of the density p. For the S model 
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we obtain 


rKiax 

( v )s — l V [(1 — P) ^min-ft;*^*} (^) + P fimin{v*+Av s ,V mEiX }{ v )] 
Jo 


(19) 


= (1 — P) min{w*,u*} + P 


v * 


Avg, 


if v* + Avg < V m 


V* “t - (Vmax U*) , if V* p Avg V m 


In contrast, if we consider the \ model we have 


(v) v = v 


/''I D\ X . f..\ i D X[w,,min{i;»+Ai> x ,V r max }]( ?; ) 

yl -P ) 0 m [ n f v v *\{V) -\~ r r A T/ d 

mm {?;* + Av x , l/ max } - u* 


dw 




= (1 — P) min-{P*, t;*} + P < 


Av 


tidr;, if v* + Av x < V max 


X Ji, 


„ V n ;■] x C* 


V dv, if U* +A'« X > Fmax 


and thus 




Ally 


if ?;* + Aiu, < V n 


(v) x = (1 — P) min-{P*, u*} + P 


X — v max 


V* P 2 (^max ^*); if P Av x > Vmax 


( 20 ) 


By comparing the last lines of (191 and (201, it is clear that 


( v ) x = (v)g VP* < V max — Au x , provided Avg = \Av x . 


( 21 ) 


Remark 11 (Uniformly distributed deceleration). The computation of the expected speed shows 
a link between the <5 model and the x model. In fact, under the constraint Avg = ^Av x , the two 
models provide the same expected speed. Similarly, if we consider a braking scenario in which 
the candidate vehicle brakes to a speed uniformly distributed in an interval centered on v* , then 
the expected speed results to be again (1 — P)v*. 

Remark 12. Let us compare the A : x matrices (Figure |9j) for r < j < N — r with a given Av and 
the corresponding A 3 § matrices (Figure [ 3 ]) with For the case r £ N, the isolated nonzero row 
of A J S is at j — |, which corresponds to the middle of the green shaded area in A 3 . Moreover, for 
any fixed Av, it can be proved that, for 5v —> 0, the sum of the quantities located in the shaded 
area of A 3 X is equal to the total contribution provided by the [j — |)-th row of the S matrices 
obtained with the acceleration parameter . In other words, as Sv —> 0, the total effect of the 
(j — §)-th row of A 3 S with ^ is spread over r + 1 rows in the matrices A 3 X with Av, which are 
the rows shaded in green in Figure [9] 
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5 Macroscopic properties 


Macroscopic acceleration In order to explain the relation of Au with the acceleration of the 
vehicles in the model, we compute the rate of change of the macroscopic velocity: 


du 

Ft 


d_ 

Ft 

V 

P Uv 

V 


~ [ vf(t,v)dv = - [ vQ[f,f](t,v)dv 

P Jv J P Jv 

[ [ A(u*->u|u*;p)/*/*du*du* - p [ u/(t,u)du 
Jv Jv Jv 

[ [ (v)f*f*dv*dv*-p[ vf(t,v)dv 

Jv Jv Jv 


P Uv Jv 

where (u) is a function of u* and u* as defined in (18I . 


du 

Ft 


In the case of the S model, ( v)s is given by (|19|) and thus 

(1 ~P) 

+ Pp I 


fV m 


f ‘ ^max Al) 


v* /*/* du*du* + (1 — P) 

rVmz 


u*/»/* du*du* 


(v* + A vg)f* du* + Pp 


hmax/* du* p 


vf(t, v) du 


'Vmax-Au 


Given an initial distribution f(t = 0,u), the equation above yields the evolution of the macro¬ 
scopic acceleration in time. It is easy to study analytically this quantity at the initial time. In 
particular, we compute the initial acceleration in the case in which all vehicles are still but the 
density is below the value for which P = 1/2 (that is p = 1/2 when taking P = 1 — p). By 
considering an initial distribution with all vehicles in the lowest velocity class, i.e. of the form 
/(0,u) = I£xii(u), we have 


du 

Ft 


= r/PAvs 


Knax Av$ 


/* dr>* + O(Sv) = rjpPAvs + O(Sv). 


t —o 


The above equation shows that the acceleration of the vehicles in the S model depends linearly 
on Au. Analogously, for the \ model, using (201, one obtains 


du 

Ft 


t=o 


= \pP p ^v x + O(Sv) 


which reinforces the remark made in (21) about the similarities of the \ an d the <5 model when 
Au x = 2 Au <5 . 

Remark 13 (Acceleration). Recall that (??p) _1 is a time. Thus rjpAv is the built-in acceleration 
of the model, which, not surprisingly, is linked to Au. We can use dimensional arguments to 
estimate the order of magnitude of Au. According to Lebacque [16] , the maximum acceleration 
of cars is approximately aLB = 2.5 m/s 2 . The maximum speed is approximately R m ax — 28m/s, 
and we expect the maximum acceleration when P = 1. Thus ppAv — aLB- 

The estimates above provide the trend of the macroscopic acceleration starting from rest. 
For the general case, we now study the evolution of the macroscopic velocity u in time, up 
to steady state. These data are shown in Fig. [lO] and Fig. m for various combinations of 
the model parameters. The results shown are obtained integrating the equations for the <5 and 
the x model found in (131 and respectively, with r € N, and computing at each time 

u(t) = lJ v vf(t,v) du. 
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Figure 10: Evolution of the macroscopic velocity in time. Left: comparison of different values of 
T and 5v. The dot-dashed lines without markers correspond to the \ model. Right: relaxation 
to steady state for different combinations of 77 and T. 


p=0.65 p=0.9 




Figure 11: Evolution of the macroscopic velocity in time, for different values of T and rj. Left: 
p = 0.65. Right: p = 0.9. 


Figure [lO] shows a typical case in which we expect acceleration. The density is p = 0.15, 
well below the value corresponding to P = 1/2 when P = 1 — p, and we start with an initial 
distribution in which /i(f = 0) = p, while fj(t = 0) = 0, for all j > 2. Thus initially all vehicles 
are still, and, since the density is low, they will accelerate to reach the maximum speed. The 
duration of the transient depends on the product r]Av = 77 /T, for a fixed density, as is apparent in 
the right panel of the figure, because the acceleration, i.e. the slope of the curves, is proportional 
to i]Av. The left panel shows the effect of the grid discretization, i.e. the role of Sv = 1/(N — 1). 
It is clear that the discretization grid has no influence on the results, as expected from Theorem 
[6] The dot-dashed lines without markers show the evolution of the macroscopic velocity obtained 
with the x model. The colour code is chosen to ensure that the curves with Av$ = \Av x are 
drawn in the same colour. As expected, the macroscopic velocity for the \ and the (5 model 
behave very similarly, provided the parameter Av is chosen correctly. 

Next, in Figure we show the evolution of the macroscopic velocity in two cases when 

we expect deceleration for the S model. Namely, we consider p = 0.65 in the left panel and 
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p = 0.9 in the right panel. The initial distribution is /jv(t = 0) = p — e, f±(t = 0) = e, and 
fj(t = 0) = 0, j = 2,... ,N — 1. The value e is introduced to ensure convergence to the stable 
equilibrium, see Remark [9] Here e = 0.01. In other words, we start with a congested traffic, 
in which initially almost all vehicles are traveling at the fastest speed available. Clearly, this 
situation is somewhat artificial, but surely we expect the vehicles to brake. Since braking does 
not depend on Av, we expect that the relaxation time towards equilibrium depends mainly on 
77 and only weakly on T. This is clearly seen in both pictures. The macroscopic speed to which 
the model relaxes on the other hand will depend on Av and on p, but not on rj. Note that when 
p = 0.9, in all cases considered here, the equilibrium speed is nearly zero: in fact the traffic is 
extremely jammed. For p = 0.65 instead, we expect that the traffic will have a residual speed, 
because we are well below the value P = 1/2, but cars are not “bumper to bumper”, and this 
residual speed does depend on Av. 


Fundamental diagrams As already discussed in the previous section, the nonzero elements of 
the matrix A 3 X can be lumped in the matrix Agj for N sufficiently large, with the only exception 
of the elements in the rectangle r x N (see Figure [9| of the matrices for j = N — r + 1,..., N. 
This is shown, for instance, in the evaluation of the expected value of the resulting speeds due 
to acceleration interactions in (19 ) and (p0| 
values close to W 


which are comparable, except again for high speed 
(and different from it by at most Av). Thus, although the x model is 


apparently more refined then the S model, we expect both models to provide similar macroscopic 
information, for large N. This is usually analyzed by computing the density and the flux as 
moments of the asymptotic kinetic distribution f°° ( v ): 


P = 


/' 

Jo 


f°°(v)dv, (pu) = / vf°°(v)dv 
/0 Jo 

and by studying the characteristics of the related fundamental diagram which is obtained plotting 
the flux against the density. 

Notice that, for the S model, Theorem [G] ensures that only few velocities, obtained with 
Sv = Av, are necessary to describe completely the exact asymptotic kinetic distribution. We 
expect therefore that the macroscopic behavior of the (5 model will be apparent even on the 
coarse velocity grids, i.e. for r = 1. 


Figure 12 shows the fundamental diagrams provided by the S model (blue curves) and the 
X model (red curve), computed with A vs = \Av x , for two different values of A vs. In the left 
panel, r = 1 , while r = 20 on the right. The figure shows that the diagram of the x model is 
very similar to the diagram of the S model when TV —> 00 and this result is in agreement with 
the fact that the expected output speed of the two models is mostly the same (i.e., the same in a 
large range of pre-interaction speeds) when choosing the acceleration parameter of the S model 
as a half of the acceleration parameter of the x model. The only difference is provided by the 
maximum speeds which, as already noted, are slightly different. Note that the similarity of the 
fundamental diagrams does not mean that the asymptotic equilibrium functions of the x and of 
the S model converge to the same function as N tends to 00. 

Observe that the fundamental diagrams given by the S model in both plots in each line of 
Figure 12 use the same information. In fact, following the results of Theorem[6j the macroscopic 
flux is given by 

N n T+l 

r v max _ 1 j 1 

Fluxj(r) = / vfff{v)dv = X^iTTj / vdv = '52( f 0i v (i-i)r +1 

jo j—i Mi I •' b 


2=1 


where Vj denotes the center of the cell 


system (131 with Av/Sv = r € N. 


f r is the vector containing the equilibria of the 
Recalling the definition of Ij, we have that v\ = Aw/4r, 


Ij and 
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Fundamental diagrams of traffic, r=1 Fundamental diagrams of traffic, r=20 




Fundamental diagrams of traffic, r=1 Fundamental diagrams of traffic, r=20 




Figure 12: Fundamental diagrams resulting from the <5 model (blue *-symbols) and from the \ 
model with acceleration parameter Av$ = | Av x (red circles). The dashed line is the flux of the 
S model in the limit r —> oo. 


i>n = Knax — Av/Ar and up_ 1 ) r+1 = {l — l)Az>. Thus, in order to compute the fundamental 
diagram of the S model with any value of r, it is enough to compute the equilibria fi using r = 1 
and then compute the flux with the formula above. In particular, using only fi, one may also 
compute the fundamental diagram of the S model also in the limit r —> oo with the formula 


T+l 


Fluxs(oo) = ^(fi)z(Z - l)Al 


1=1 


The dashed blue line in all panels of Figure 12 shows the quantity Fluxa(oo) just defined. Note 


that in the case of the x model, for each value of r, one has instead to compute the full equilibrium 
distribution with N = rT + 1 velocities. 

When increasing r, we observe that the flux at p max approaches zero. This is because for 
p = Pmax) (fi)i is the only non zero component at equilibrium, all vehicles travel at a velocity in 
the lowest speed class I\ and the flux is therefore ^(fi)i. Similarly, in the free phase all vehicles 
travel at a velocity in the highest speed class In and the flux is therefore (Knax — : fy)(fi)T+i = 
(Knax — ^)p- The free-phase flux is therefore linear in p and its slope approaches Vj na x when 
r —► oo. 

In Figure |i~2} we observe that both models provide a sharp decrease in the flux, beyond the 
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8 model - Fundamental diagrams of traffic, r=1, T=4 8 model - Fundamental diagrams of traffic, r=20, T=4 




Figure 13: Fundamental diagrams resulting from the S model with acceleration parameter A vg = 
The probability P is taken as in ([5]) with 7 = 1 (blue data), 7 = 3/4 (green data) and 7 = 1/4 
(cyan data). The dashed lines are the fluxes in the limit r —> 00 . 


critical density , namely the value of the density marking the transition from free to congested 
flow. This phenomenon is well known in traffic modeling, and it is called capacity drop , see |28] 
and references therein. From Theorem [ 6 ] it is apparent that, for the S model, the critical density 
corresponds to a bifurcation of the equilibrium solutions. In fact, one deduces that for P > | the 
equilibrium distribution is f°°(v) = p<V max (rO, which means that all vehicles travel at maximum 
speed. Only when P < \ the lower speed classes begin to fill up. Thus, the physical concept 

of phase transition in traffic flow theory has a rigorous mathematical counterpart in the present 

JL 

model. Using the law given in ^ the value of p for which P = 1/2 is p c := (Z ) 7 and then we 
may act on 7 in order to change the critical density. For instance, see Figure [13] in which we plot 
three fundamental diagrams of the 5 model with P as in ([ 5 ]) and 7=1 (^-markers), 7 = 3/4 
(x), 7 = 1/4 (A). 


Remark 14. Observe from Figure [12] that the critical density of the \ model approaches the 
critical density of the S model when N —> 00 . In fact, since the matrix A* —> A\ for r 00 , we 

(*?)r 


also have that (f£)j 


A$ for r - 

More precisely, the analogous of (|3.1| for the x model is 


- (1 - P) fl + 


(l- 2P+£)„,.=<, 


and the stable equilibrium is thus 


0 

p^ + oO 


P > 


2 r 

4r—1 


otherwise 


In Figure[l4]we show the fundamental diagrams of the x model for r = 1 and r = 20, together 
with a few representative /./s at equilibrium, as functions of p. In the left part, for r = 1, two 
phase transitions appear in the fundamental diagram (top left). Comparing with the bottom 
left plot, the origin of this phenomenon can be appreciated. A first transition occurs when the 
density becomes large enough to force a few drivers to brake: thus the second largest speed class 
In-i starts being populated (green dashed curve), while the fastest speed class begins to be 
depleted (red curve). A second transition occurs when some vehicles enter the lowest speed class 


27 














r=1 



Density p 


r=20 



Figure 14: Top: fundamental diagrams provided by the \ model with N = 4 (left) and N = 61 
(right) velocities. Bottom: equilibria of the function /i (blue solid line), /n-i (green dashed) 
and /at (red dot-dashed) for any density in [0,1]. 


(blue curve). This latter transition is the one that, when increasing r, moves towards the critical 
density p = 1/2, see Remark 14 The first phase transition is not observable for large r, because 
fN—i is related to the velocity vn-i -* Vmax, as Sv —> 0, so that the transition of vehicles from 
In to In —i is not enough to determine an abrupt change in the flux. 


Comparison with experimental data Figure 15 shows the comparisons of the results pro¬ 
duced by the S model with experimental data published in m In the left figure we have tuned 
the critical density to reproduce the correct position of the phase transition. The experimental 
data are normalized and the fundamental diagram computed by the model is provided for all 
values of the density between 0 and p max , which corresponds to a situation in which all vehi¬ 
cles are bumper-to-bumper and still. The figure on the right stems from the observation that 
experimental data contain a residual movement even in the congested phase. Thus the bumper- 
to-bumper situation is never actually observed. Therefore in the figure on the right we also 
tune the maximum density p max actually observed, with p max < p ma x- In this case we obtain a 
very good agreement with experimental data. With the present model we do not reproduce the 
scattering of the data, but this can be explained keeping into account a mixture of two different 
populations of drivers and/or vehicles, as we have proposed in [251 . 


6 Conclusions and perspectives 

In this work we have studied two kinetic models for vehicular traffic based on a Boltzmann-like 
term describing binary interactions. We have assumed a continuous space of microscopic speeds 
and we have analyzed the space homogeneous case to study the asymptotic behavior of the 
distribution function together with the resulting flux-density diagrams. 
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Figure 15: Comparison between experimental data and the diagram resulting from the S model, 
with An = 1/4, P = 1 — p 1 / 4 . The experimental diagram is reproduced by kind permission of 
Seibold et al. m- 


Our models are characterized by a parameter Aw, that has physical relevance and is related 
to the maximum speed variation in a unit of time. The models are defined by the transition 
probability of gaining a given velocity and they differ only in the modeling of the acceleration 
interaction. 

First of all we have studied the case in which the resulting speed after an acceleration is 
obtained by a velocity jump from w* to w* + Av, where w* is the pre-interaction speed. We 
have referred to this model as S model and we have found a class of asymptotic distributions 
which is atomic with respect to the velocity variable. In other words it is a linear combination of 
Dirac delta functions centered in a finite number of velocities. The number T of delta functions 
contributing to the stable equilibrium distribution is controlled by the acceleration parameter 
through the relation Av = V ma , X /T. This result means that the number of discrete velocities 
necessary to completely describe the equilibrium distribution function is implicitly determined 
by the acceleration parameter Av and therefore is small. 

Instead, in the \ model we have prescribed the acceleration interaction in a way that is closer 
to the modeling given in m, but again respecting the physical relevance of Av. In fact, we 
have assumed that the output speed after acceleration is uniformly distributed over the range 
[w*, w* + Aw]. We have shown that the \ model with acceleration parameter Aw gives a macro¬ 
scopic behavior similar to the one provided by the simpler S model with acceleration parameter 
Aw/2, as it can be seen by studying the macroscopic properties of the two models and comparing 
their fundamental diagrams, notwithstanding the fact that the respective asymptotic distribution 
functions do not approach each other. Thus the \ model, despite the more sophisticated descrip¬ 
tion of interactions, gives the same macroscopic information of the simpler and computationally 
much cheaper S model, at least at equilibrium. Moreover we have proved that both models 
provide a bounded macroscopic acceleration, studying the evolution in time of the macroscopic 
velocity, and its relation with the parameter Aw. 

The results obtained in this work suggest that a small number of velocities is sufficient for 
the kinetic modeling of traffic. This is crucial to make kinetic modeling of complex traffic 
flows amenable to computations. Note also that here the acceleration remains controlled by 
Aw, in contrast to models based on a lattice of velocities, in which the possible outcomes of an 
interaction and the acceleration of vehicles depend on the particular lattice chosen [B]. Moreover 
the complete knowledge of the equilibrium distribution is crucial to derive macroscopic models 
with a rich enough closure law resulting naturally as consequence of the microscopic interactions. 
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Thus, without the need of prescribing heuristic speed-density relations, we obtain fundamental 
diagrams with a phase transition and a capacity drop as those occurring in experimental data. 

Finally, the particular structure of the equilibrium distribution provided by the <5 model allows 
one to generalize this framework to the case of multipopulation models, as in (24j . in order to 
recover multivalued fundamental diagrams, see also [25] . 
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A Matrix elements for the discretization of the x model 

In order to compute the Af matrices, we observe that the gain operator of the S model given 
in (|9]) differs from the gain operator of the \ model only in the terms proportional to P appearing 
in the equation (16). Therefore, we just show the terms resulting from 




rVmax-All 


f-Vm 


X[v„v,+Av](v)f*dv* + Pp 


I'max Av 


ax] (' t 0 

Vmax C* 


f*dv*. 
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When the terms above are integrated over the cells 1 1 , we get 


[ -G[f,f](t,v)dv=^h. 

Jh V 4r 


For j = 2,..., r, 

[ 1 G[f,f](t,v)dv= F ^Y / f h + ^f j . 

•M r ti 2r 

For j = r + 1, 






h-j—r-\-l 


For j = r + 2,..., N — r — 1, 


^G[/, /](t,u)dv=^/j_ r 



h=j— r+1 


For j = N — r, 



G[/, /](*, w)dv =^-fj_ r 


Pp 


j -1 

E A + p P 

h=j-r-\-l 


3 

8 r 


1 

2 


- — r ) !og 


2 r 


2r- 


For j = N — r + 1,..., N — 1, 

t-G\s,m^ = ^rfi-r+— E 

It. 7 ? 2 r r ' 

^ J-n I U — 


AT— r— 1 


h=j-r -\-1 


1 , / 2r 

2? + l0g (s-I 


fN-r 


4-1 /TV — ft + | \ 

+ +J+ log (^A) A 

+ P,,[l + (i + i-iv)lo g (+j±l 




Finally, for j = N, 

f 1 


//iv n 


~~G[f, f](t,v)dv = Pp 


1 !, / 2r 

8 r + 2 bg ( 2r - 1 


fj—r + 


j -1 

E lo § 


Pp +1 ,_,^-ft+i 


h=j— r+1 


N-h-\ 


The matrices A J can be formed by removing the underlined terms in (121 with r 


adding the contributions given in (A.l). 


(A.la) 


(A. lb) 


(A.lc) 


(A.ld) 



(A.le) 

(A.lf) 


h + PpfN- 
(A.lg) 

£ N and 
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